#Nature's Kidneys: the Role of Wetland Reserve Easements in Restoring Water Quality
#Nicole Karwowski and Marin Skidmore
#Replication File
#2025

#Clean and prep the geospatial NRCS shapefiles

#open libraries
library(ggplot2)
library(raster)
library(sf)
library(dplyr)
library(readr)
library(mapview)

##########################################################################################

#Clean WRP shapefile data
#open wetland easements sf
setwd("C:/Users/16308/OneDrive - Montana State University/UW_Madison/Water quality/Replication Package")

ease<-st_read("Data/NRCS easements/Raw/NRCS_easements_sf/easements.shp")

colnames(ease)
table(ease$NEST_PROGR)

#KEEP
#ACEP-WRE: Agricultural Conservation Easement Program- Wetland Reserve Easement
#EWPP-FPE: Emergency Waterhsed Propetection Program- Floodplain Easement
#RCPP-WRE: Regional Conservation Partnership Program-Wetland Reserve Easement
#WRP: Wetland Reserve Program

#DELETE
#ACEP-ALE: Agricultural Conservation Easement Program-Agricultural Land Easement
#FRPP: Farm and Ranch Protection Program
#GRP: Grassland Reserve Program 
#HFRP: Healthy Forest Reserve Program
#OSL: Other Stewardship Lands
#RCPP-ALE:  Regional Conservation Partnership Program- Agricultural Land Easement
#RCPP-HFRP: Regional Conservation Partnership Program- Healthy Forest Reserve Program

#keep ACEP-WRE, EWPP-FPE, EWRP, RCPP-WRE, WRP. 
ease<-subset(ease, ease$NEST_PROGR=="ACEP-WRE" | ease$NEST_PROGR=="EWPP-FPE" |
               ease$NEST_PROGR=="EWRP" | ease$NEST_PROGR=="RCPP-WRE" | 
               ease$NEST_PROGR=="WRP")

table(ease$NEST_STATE)
table(ease$Admin_Area)
table(ease$Method)
table(ease$NEST_EASEM)

#make sf valid
ease<-st_make_valid(ease)

list(ease)
#create year variable
ease$year<-substring(ease$NEST_C_CLO, 1, 4)

#create month var
ease$month<-substring(ease$NEST_C_CLO, 6, 7)

#create year-month variable
ease$yearmonth<-paste0(ease$year, ease$month, sep="")

ease<- ease %>% relocate(year, .before=NEST_C_CLO)
ease<- ease %>% relocate(month, .after=year)
ease<- ease %>% relocate(yearmonth, .after=month)

#replace NEST_EASEM
ease$NEST_EASEM<-ifelse(ease$NEST_EASEM=="30-Year", "30 Year Easement", ease$NEST_EASEM)
ease$NEST_EASEM<-ifelse(ease$NEST_EASEM=="Permanent", "Permanent Easement", ease$NEST_EASEM)

#create binary for permanent
ease$permanent<-ifelse(ease$NEST_EASEM=="Permanent Easement", 1, 0)
ease<- ease %>% relocate(permanent, .after=NEST_EASEM)

#look into the timing
table(ease$year)
table(ease$month)

#create unique easement id
ease<-ease[order(ease$NEST_C_CLO), ]
ease<-mutate(ease, EasementID=row_number())
ease<- ease %>% relocate(EasementID, .before=OBJECTID)

#delete some vars and rename some vars
names(ease)
ease<-ease[c("EasementID", "OBJECTID", "Admin_Area", "CalcAcres", 
             "NEST_STATE", "NEST_Count", "NEST_PROGR","permanent", "year",
             "month", "yearmonth", "NEST_C_CLO", "Shape_Area", "geometry")]

ease<- ease%>% relocate(CalcAcres, .before=Shape_Area)

names(ease)<-c("EaseId","NrcsId","AdminArea", "State", "County", "Program",
                "Permanent", "Year", "Month", "YearMonth", "Closing", 
               "CalcAcres", "ShapeArea", "geometry")
names(ease)
list(ease)

#keep permanent easements only
ease<-subset(ease, ease$Permanent==1)


######################################################################

#merge it with the easement detailed data to get the restoration date
load("Data/NRCS easements/Processed/EasementsDetailedCleaned.RData")

#turn detailed data into sf
Details<-st_as_sf(Details, coords=c("Longitude", "Latitude"), crs=crs(ease))
Details<-subset(Details, select=c(Id, Program, State, CountyName, FinalAcres, ClosingDate, ClosingYM, RestorationDoneYM, geometry))
crs(ease)==crs(Details)

#Merge lat and long with geometry
merge<-st_join(ease, Details, join=st_contains, left=FALSE)
length(unique(merge$EaseId))

#keep matched obs and only distinct ones
good_merge<-subset(merge, !is.na(merge$Id) & !duplicated(merge$EaseId))

########################################################################
#Merging in the second way

#collect the leftovers
redo<-subset(ease, !(ease$EaseId %in% good_merge$EaseId))
unmatched<-subset(Details, !(Details$Id) %in% good_merge$Id)

#try to do merge with closing date, county name
redo$Closing<-substring(redo$Closing, 1, 10)
redo$Closing<-as.Date(redo$Closing)

#and calculated acreage
redo$Acres<-as.integer(redo$CalcAcres)
unmatched$Acres<-as.integer(unmatched$FinalAcres)

#trim the county names
redo$County<-str_trim(redo$County)
unmatched$CountyName<-str_trim(unmatched$CountyName)


merge2<-merge(st_drop_geometry(redo), st_drop_geometry(unmatched), 
              by.x=c("Closing", "County", "Acres"), 
              by.y=c("ClosingDate", "CountyName", "Acres"),
              all.x=FALSE, all.y=FALSE)


length(unique(merge2$EaseId))

#only keep the nonduplicated matches
good_merge2<-subset(merge2, !duplicated(merge2$EaseId) & !duplicated(merge2$Id))


###############################################################################
#Merge #3
#again, compile the leftovers. we have 2,103 left that we want to match
redo1<-subset(redo, !(redo$EaseId %in% good_merge2$EaseId) )
unmatched1<-subset(unmatched, !(unmatched$Id) %in% good_merge2$Id)

unmatched1$Year<-substring(unmatched1$ClosingYM,1,4)

#match on year, county, acres
merge3<-merge(st_drop_geometry(redo1), st_drop_geometry(unmatched1), 
              by.x=c("Year", "County", "Acres"), 
              by.y=c("Year", "CountyName", "Acres"),
              all.x=FALSE, all.y=FALSE)

length(unique(merge3$EaseId))

#only keep the nonduplicated matches
good_merge3<-subset(merge3, !duplicated(merge3$EaseId) & !duplicated(merge3$Id))


###############################################################################
#Another attempt to merge
redo2<-subset(redo1, !(redo1$EaseId %in% good_merge3$EaseId) )
unmatched2<-subset(unmatched1, !(unmatched1$Id) %in% good_merge3$Id)

#try rounding rather than the integer
redo2$Acres<-round(redo2$CalcAcres)
unmatched2$Acres<-round(unmatched2$FinalAcres)

#match on year, county, acres
merge4<-merge(st_drop_geometry(redo2), st_drop_geometry(unmatched2), 
              by.x=c("Year", "County", "Acres"), 
              by.y=c("Year", "CountyName", "Acres"),
              all.x=FALSE, all.y=FALSE)

length(unique(merge4$EaseId))

#only keep the nonduplicated matches
good_merge4<-subset(merge4, !duplicated(merge4$EaseId) & !duplicated(merge4$Id))


###############################################################################
#Another attempt to merge
redo3<-subset(redo2, !(redo2$EaseId %in% good_merge4$EaseId) )
unmatched3<-subset(unmatched2, !(unmatched2$Id) %in% good_merge4$Id)

#state and closing date?
merge5<-merge(st_drop_geometry(redo3), st_drop_geometry(unmatched3), 
              by.x=c("YearMonth", "County"), 
              by.y=c("ClosingYM", "CountyName"),
              all.x=FALSE, all.y=FALSE)
merge5$diff<-abs(merge5$CalcAcres-merge5$FinalAcres)
summary(merge5$diff)

#only keep the nonduplicated matches
good_merge5<-subset(merge5, !duplicated(merge5$EaseId) & !duplicated(merge5$Id) & merge5$diff<5)

###############################################################################
#Another attempt to match
redo4<-subset(redo3, !(redo3$EaseId %in% good_merge5$EaseId) )
unmatched4<-subset(unmatched3, !(unmatched3$Id) %in% good_merge5$Id)

#nearest neighbor and closing date month match
nn<-st_join(redo4, unmatched4, join=st_nearest_feature, left=FALSE, maxdist=500)

nn1<-subset(nn, nn$Year.x==nn$Year.y)
nn1<-subset(nn1, nn1$State.x==State.y)
nn1<-subset(nn1, !is.na(nn1$RestorationDoneYM))
nn1<-subset(nn1, nn1$County==nn1$CountyName)

#acreage differential
nn1$diff<-abs(nn1$CalcAcres-nn1$FinalAcres)
summary(nn1$diff)

nn1<-subset(nn1, nn1$diff<10)
length(unique(nn1$EaseId))

#try year and acres
good_merge6<-nn1
good_merge6<-subset(good_merge6, select=c(EaseId, Id, RestorationDoneYM))
good_merge6<-st_drop_geometry(good_merge6)


################################################################################
redo5<-subset(redo4, !(redo4$EaseId %in% good_merge6$EaseId) )
unmatched5<-subset(unmatched4, !(unmatched4$Id) %in% good_merge6$Id)

#county and year
merge7<-merge(st_drop_geometry(redo5), st_drop_geometry(unmatched5), 
              by.x=c("Year", "County"), 
              by.y=c("Year", "CountyName"),
              all.x=FALSE, all.y=FALSE)

#only keep the nonduplicated matches
good_merge7<-subset(merge7, !duplicated(merge7$EaseId) & !duplicated(merge7$Id))
good_merge7<-subset(good_merge7, good_merge7$State.x==good_merge7$State.y)

#################################################################################
#create a key
key<-subset(st_drop_geometry(good_merge), select=c(EaseId, Id, RestorationDoneYM))

key<-rbind(key, subset(good_merge2, select=c(EaseId, Id, RestorationDoneYM)))

key<-rbind(key, subset(good_merge3, select=c(EaseId, Id, RestorationDoneYM)))

key<-rbind(key, subset(good_merge4, select=c(EaseId, Id, RestorationDoneYM)))

key<-rbind(key, good_merge6)

key<-rbind(key, subset(good_merge7, select=c(EaseId, Id, RestorationDoneYM)))

#97% we were able to match

length(unique(key$EaseId))

#if unable to get the merge, assume the restoration date
##on average, it takes 913 days between closing and restoration. 
#that is equal to 30 months, or 2.5 years. 

ease1<-merge(ease, key, by=c("EaseId"), all.x=TRUE, all.y=FALSE)

ease1$ClosingDate<-as.Date(ease1$Closing)
ease1$RestorationDate<-ease1$ClosingDate+913
ease1$RestorationDate<-as.character(ease1$RestorationDate)
ease1$RestorationDate<-substring(ease1$RestorationDate, 1, 7)
ease1$RestorationDate<-gsub("-", "", ease1$RestorationDate)

ease1$RestorationDoneYM<-ifelse(!is.na(ease1$RestorationDoneYM), ease1$RestorationDoneYM,
                                ease1$RestorationDate)

table(ease1$RestorationDoneYM==ease1$RestorationDate)

#save it!
ease1<-subset(ease1, select=-c(ClosingDate, RestorationDate))
colnames(ease1)

st_write(ease1, "Data/NRCS easements/Processed/NRCS_easements_sf_clean/NRCS_easements_sf_clean.shp")
